# 1. Packages qui interfacent des données OSM
# install.packages('maposm', repos = 'https://riatelab.r-universe.dev')
library(tidygeocoder) # géocodage
library(maptiles) # import de tuiles OSM (raster)
library(osmdata) # import de données OSM (vecteur)
library(maposm) # import de données OSM (couches géo)
library(osrm) # calcul d'itinéraires
# 2. Autres packages utilitaires
library(sf) # manipulation de données vectorielles
library(terra) # manipulation de données raster
library(mapsf) # cartographie thématique
library(maplegend) # légendes
library(spatstat) # analyse de semis de points
library(stplanr) # segmentiser des lignes
library(mapview) # cartographie interactive
library(lwgeom)Utiliser OpenStreetMap avec R - mise en pratique
0. Packages
1. Géocodage
Le package tidygeocoder (Cambon et al. 2025) permet d’utiliser des services de géocodage en ligne. Par défaut, le package utilise Nominatim qui se base sur les données OSM.
pt <- data.frame(address = "28 Place de l'Hôtel de Ville, 60200 Compiègne")
pt <- geocode(.tbl = pt, address = "address", quiet = TRUE) |>
st_as_sf(coords = c("long", "lat"), crs = 4326)
ptSimple feature collection with 1 feature and 1 field
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 2.826361 ymin: 49.41777 xmax: 2.826361 ymax: 49.41777
Geodetic CRS: WGS 84
# A tibble: 1 × 2
address geometry
* <chr> <POINT [°]>
1 28 Place de l'Hôtel de Ville, 60200 Compiègne (2.826361 49.41777)
2. Import de tuiles OSM
Le package maptiles (Giraud 2025) permet de télécharger des tuiles chez différents fournisseurs, dont OpenStreetMap, au format raster. 20 niveaux de zoom sont disponibles, selon la taille de l’objet géographique à cartographier. Un zoom de niveau 10 correspond à une échelle 1/500.000.
Plusieurs modèles de tuiles sont disponibles ici, certains nécessitent une inscription.
emprise <- st_buffer(pt, 3000) |>
st_bbox()
tiles <- get_tiles(x = emprise,
project = FALSE,
crop = TRUE,
cachedir = "cache",
zoom = 14)
pt_reproj <- st_transform(pt, crs = "EPSG:3857")
mf_raster(tiles, expandBB = c(rep(-.2,4)))
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_annotation(pt_reproj, txt = "Mairie de Compiègne")
mf_title("Tuiles OSM - Défaut")# Autre exemple de fournisseur
tiles <- get_tiles(x = emprise,
project = FALSE,
crop = TRUE,
provider = "OpenStreetMap.HOT",
zoom = 14)
mf_raster(tiles, expandBB = c(rep(-.2,4)))
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_annotation(pt_reproj, txt = "Mairie de Compiègne")
mf_title("Tuiles OSM - OpenStreetmap.HOT")3. Import de couches géographiques OSM
Pour la cartographie, le package maposm (pas sur le CRAN). Il se base sur les fonctions du package osmdata (Padgham et al. 2023) (partie suivante) pour télécharger les données OSM, au format vectoriel.
La fonction principale, om_get(), permet de récupérer certaines couches géographiques utiles (routes, parcs, bâti, etc), autour d’un point, dans un rayon donné.
res <- om_get(x = c(st_coordinates(pt)[1],
st_coordinates(pt)[2]),
r = 2000)Getting urban areas: 0.84 sec elapsed
Getting buildings: 17.2 sec elapsed
Getting green areas: 1.58 sec elapsed
Getting roads: 0.47 sec elapsed
Getting streets: 1.33 sec elapsed
Getting railways: 0.48 sec elapsed
Getting water bodies: 2.78 sec elapsed
names(res)[1] "zone" "urban" "building" "green" "road" "street" "railway"
[8] "water"
mf_map(res$zone, col = "#f2efe9", border = NA, add = FALSE)
mf_map(res$urban, col = "#e0dfdf", border = "#e0dfdf", lwd = .5, add = TRUE)
mf_map(res$green, col = "#c8facc", border = "#c8facc", lwd = .5, add = TRUE)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$building, col = "#d9d0c9", border = "#c6bab1", lwd = .5, add = TRUE)
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("Autour de la mairie")
mf_arrow()
mf_scale(size = 500, scale_units = "m")
mf_credits("OpenStreetMap contributors, 2025", pos = "bottomleft")5. Carroyage
Une fois les points d’intérêt récupérés, il est possible de les agréger dans une maille.
grid <- st_make_grid(res$zone, cellsize = 300, square = FALSE)
grid <- st_intersection(grid, res$zone)
grid <- st_sf(ID = 1:length(grid), geom = grid)
# compter
inter <- st_intersects(grid, poi_reproj, sparse = TRUE)
grid$n_poi <- lengths(inter)
om_map(res, title = "Autour de la mairie")
mf_map(grid, var = "n_poi", type = "choro", breaks = c(0,1,3,5,10,15),
alpha = .6, border = NA, leg_title = "Nombre \nd'aménités",
leg_val_rnd = 0, leg_pos = "topright", add = TRUE)
mf_title("Aménités gustatives - données carroyées")6. Lissage (KDE)
Le package spatstat (Baddeley, Turner, and Rubak 2025) regroupe de nombreuses fonctions pour analyser des semis de points. Il utilise (génère ?) des objets de classe ppp (planar point pattern), dans des fenêtres d’observations, de classe owin (observation window).
La fonction density.ppp() calcule des valeurs à partir d’une estimation par noyau (KDE). Elle prend en entrée une portée (sigma), ainsi qu’une résolution au sein de laquelle sera calculée la densité, en pixels (eps).
p <- as.ppp(X = st_coordinates(poi_reproj), W = as.owin(res$zone))
ds <- density.ppp(x = p, sigma = 100, eps = 10, positive = TRUE)
# Calcul densité de restaurants par hectare
r <- rast(ds) * 100 * 100
crs(r) <- st_crs(poi_reproj)$wkt
rclass : SpatRaster
dimensions : 400, 400, 1 (nrow, ncol, nlyr)
resolution : 10, 10 (x, y)
extent : 312629.1, 316629.1, 6344049, 6348049 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / Pseudo-Mercator (EPSG:3857)
source(s) : memory
name : lyr.1
min value : 2.225074e-304
max value : 1.970817e+00
mf_map(res$zone, col = "#f2efe9", border = NA)
mf_raster(r, leg_title = "N. restaurants / ha.", leg_pos = "left",
pal = c( "#75871B", "#BAB97D", "#FAF3CE", "#FDE16A", "#F9B40E",
"#E88C23", "#A25933"), add = TRUE)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(poi_reproj, pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("La montagne alimentaire - Compiègne")Il est possible de convertir cette couche en format vecteur. Par défaut, la fonction as.polygons() arrondit les valeurs des pixels à l’entier.
# par défaut
r_poly <- as.polygons(r) |>
st_as_sf()
# tous les dixièmes
r_poly <- as.polygons(r, round = TRUE, digits = 1) |>
st_as_sf()
mf_map(res$zone, col = "#f2efe9", border = NA)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(r_poly, var = "lyr.1", breaks = quantile(r_poly$lyr.1), border = NA,
lwd = .2, type = "choro", leg_title = "Densité\n (n/ha)", alpha = 0.6,
pal = "Mint", add = TRUE)
mf_map(poi_reproj, pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)Par aménité, discrétisation commune
# Cafés et bars ensemble
poi_reproj$amenity[poi_reproj$amenity %in% c("bar", "cafe")] <- "bar-cafe"
p <- list()
ds <- list()
r <- list()
sel <- levels(as.factor(poi_reproj$amenity))
for(i in 1:length(sel)){
p[[i]] <- as.ppp(
X = st_coordinates(poi_reproj[poi_reproj$amenity == sel[i] ,]),
W = as.owin(res$zone))
ds[[i]] <- density.ppp(x = p[[i]], sigma = 100, eps = 10, positive = TRUE)
r[[i]] <- rast(ds[[i]]) * 100 * 100
crs(r[[i]]) <- st_crs(poi_reproj)$wkt
r[[i]] <- as.polygons(r[[i]], round = TRUE, digits = 1) |>
st_as_sf()
}
# Récupérer l'amplitude
amp <- do.call(rbind, r)
amp <- sort(c(unique(amp$lyr.1)))
bks <- quantile(amp)
pal <- mf_get_pal(n = length(bks) -1, palette = "Reds", alpha = .6, rev = TRUE)
mf_map(res$zone, col = "#f2efe9", border = NA)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(r[[1]], var = "lyr.1", breaks = bks, type = "choro", border = NA,
leg_title = "Densité\n (n/ha)", pal = pal, add = TRUE)
mf_map(poi_reproj[poi_reproj$amenity == "bar-cafe" ,],
pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("Densité de bars et cafés")mf_map(res$zone, col = "#f2efe9", border = NA)
mf_map(res$water, col = "#aad3df", border = "#aad3df", lwd = .5, add = TRUE)
mf_map(res$railway, col = "grey50", lty = 2, lwd = .2, add = TRUE)
mf_map(res$road, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(res$street, col = "white", border = "white", lwd = .5, add = TRUE)
mf_map(r[[2]], var = "lyr.1", type = "choro", border = NA, breaks = bks,
leg_title = "Densité\n (n/ha)", pal = pal, add = TRUE)
mf_map(poi_reproj[poi_reproj$amenity == "restaurant" ,],
pch = 20, cex = .2, col = "black", add = TRUE)
mf_map(res$zone, col = NA, border = "#c6bab1", lwd = 4, add = TRUE)
mf_title("Densité de restaurants")Effet de bord
7. Temps d’accès
Le package osrm (Giraud 2024) interface l’engin de routage OSRM et propose différentes fonctions pour le calcul de temps de trajet et d’itinéraires. D’après la documentation de l’API OSRM, le nombre de requêtes est limité à 512 par connection à l’API, espacées d’au moins 5 secondes.
Pour des calculs plus volumineux (échelle nationale voire européenne), il est possible d’installer une instance OSRM sur son propre serveur (voir ici).
7.1. Depuis la mairie
Pour calculer les temps de trajet (en minutes) et les distances (en mètres) entre la mairie de Compiègne et l’ensemble des restaurants, bars et cafés situés dans un rayon de 2000m, nous utilisons la fonction osrmTable(), qui présente en sortie une liste composée de deux dataframes correspondant aux coordonnées des points d’origine et de destination, et d’un vecteur correspondant à la métrique calculée. Il est possible de calculer les distances et temps de trajet entre plusieurs points en une seule requête.
# Création de la matrice O/D
o <- data.frame(X = st_coordinates(pt)[, 1],
Y = st_coordinates(pt)[, 2])
row.names(o) <- 1:nrow(o)
d <- data.frame(X = st_coordinates(poi)[, 1],
Y = st_coordinates(poi)[, 2])
head(o) X Y
1 2.826361 49.41777
head(d) X Y
1 2.820433 49.41496
2 2.824740 49.41771
3 2.824230 49.41931
4 2.823624 49.42116
5 2.826944 49.41643
6 2.832764 49.41088
# 1 requête : 1 point d'origine, 92 points de destination
mat_dist <- osrmTable(src = o,
dst = d,
measure = "distance",
osrm.profile = 'foot')
Sys.sleep(5)
mat_dur <- osrmTable(src = o,
dst = d,
measure = "duration",
osrm.profile = 'foot')
dist <- t(mat_dist$distances)
dur <- t(mat_dur$durations)
poi$dist <- dist
poi$dur <- dur
apply(poi$dur, 2, summary) 1
Min. 0.200000
1st Qu. 2.750000
Median 3.950000
Mean 6.087805
3rd Qu. 8.050000
Max. 22.700000
apply(poi$dist, 2, summary) 1
Min. 16.0000
1st Qu. 204.0000
Median 295.0000
Mean 453.6098
3rd Qu. 596.0000
Max. 1693.0000
Il est possible d’atteindre 41 restaurants, cafés ou bars en moins de 4 min à pieds depuis la mairie de Compiègne, et les trois quarts se situent à moins de 596m à pieds.
7.2. Depuis la grille
nrow(grid) * nrow(d)[1] 15498
Afin de construire des indicateurs d’accessibilité, il est nécessaire de calculer les temps de trajet vers tous les restaurants / bars / cafés de la ville depuis chaque point de grille.
Nous avons en tout 15498 couples origine / destination. Pour éviter de surcharger le serveur de démo, nous découpons nos requêtes en autant de destinations. Ainsi, chaque requête calculera le temps de trajet entre nos 189 points d’origine et un point de destination.
Les calculs ne seront pas présentés ici, et ont été effectués en local. Les résultats sont enregistrés dans le dossier data (fichier mat.rds).
Pour que l’engin de routage calcule un itinéraire, les points d’origine et de destination doivent être localisés sur une route. Si le point d’entrée ne se situe par sur une route, il sera déplacé sur le point le plus proche situé sur une route, en distance euclidienne. Cela peut avoir un impact sur l’itinéraire final ainsi que les temps de trajet, plus ou moins important selon le type d’espace (urbain ou rural).
# Coordonnées des points d'origine : centroïdes des cellules de la grille
o <- st_centroid(grid) |>
st_transform(crs = 4326)
o$X <- st_coordinates(o)[, 1]
o$Y <- st_coordinates(o)[, 2]
o <- st_set_geometry(o[, c("X", "Y")], NULL)
# Les points de destination restent les mêmes
mat <- list()
for(i in 1:nrow(d)){
mat[[i]] <- osrmTable(src = o[, c("X", "Y")],
dst = d[i , c("X", "Y")],
measure = "duration",
osrm.profile = 'foot')
Sys.sleep(5)
mat[[i]]$ID_ORIG <- c(1:nrow(o))
mat[[i]]$ID_DEST <- rep(i, nrow(o))
mat[[i]] <- data.frame("ID_ORIG" = mat[[i]]$ID_ORIG,
"ID_DEST" = mat[[i]]$ID_DEST,
"duration" = mat[[i]]$durations)
colnames(mat[[i]]) <- c("ID_ORIG", "ID_DEST", "duration")
}
mat <- do.call(rbind, mat)
head(mat)
saveRDS(mat, "data/mat.rds")Nous pouvons maintenant calculer plusieurs indicateurs d’accessibilité : temps de trajet au restaurant le plus proche, au deuxième, au troisième (question du choix), nombre de restaurants à moins de 5 minutes, 10 minutes, 15 minutes, etc.
mat <- readRDS("data/mat.rds")
# Restaurant le plus proche pour chaque point de grille
poi_min <- aggregate(duration ~ ID_ORIG, mat, FUN = min, na.rm = TRUE)
# Récupérer l'identifiant de destination associé
poi_min <- do.call(rbind,
lapply(split(mat, mat$ID_ORIG),
function(x) x[which.min(x$duration), ]))
# Cartographier
grid <- merge(grid, poi_min,
by.x = "ID", by.y = "ID_ORIG", all.x = TRUE)
om_map(res, title = "Temps d'accès au restaurant le plus proche")
mf_map(grid, var = "duration", type = "choro", breaks = seq(0,15,3),
alpha = .6, border = NA, leg_title = "Temps\n(min.)",
leg_val_rnd = 0, pal = "Greens", leg_pos = "topright", add = TRUE)# Nombre de poi à moins de 5 minutes
mat5 <- mat[mat$duration < 5 ,]
mat5 <- aggregate(duration ~ ID_ORIG, mat5, FUN = length)
grid <- merge(grid[, "ID"], mat5,
by.x = "ID", by.y = "ID_ORIG", all.x = TRUE)
grid$duration[is.na(grid$duration)] <- 0
om_map(res, title = "Restaurants à moins de 5 minutes")
mf_map(grid, var = "duration", type = "choro", breaks = c(0,1,5,10,30,40),
alpha = .6, border = NA, leg_title = "Nb de\nrestaurants",
leg_val_rnd = 0, pal = "Reds", leg_pos = "topright", add = TRUE)8. Itinéraires
Pour aller plus loin, il est possible de récupérer les tracés des itinéraires, à l’aide de la fonction osrmRoute(). Attention, l’API route d’OSRM permet de calculer des itinéraires seulement entre deux points. Une requête correspond donc à un couple origine / destination.
routes <- list()
o <- data.frame(X = st_coordinates(pt)[, 1],
Y = st_coordinates(pt)[, 2])
row.names(o) <- 1
for(i in 1:nrow(d)){
routes[[i]] <- osrmRoute(src = o,
dst = d[i,],
overview = "full",
osrm.profile = 'foot')
Sys.sleep(5)
}
routes <- do.call(rbind, routes)
saveRDS(routes, "data/routes.rds")Le package stplanr (Lovelace, Ellison, and Morgan 2024) permet initialement de modéliser et de travailler sur des systèmes de transport. Nous l’utilisons ici à des fins cartographiques, pour découper des lignes en segments.
routes <- readRDS("data/routes.rds")
routes <- st_transform(routes, crs = "EPSG:3857")
routes$n <- 1
routes$distance <- routes$distance * 1000
# La fonction overline permet de compter le nombre d'occurence des tronçons
seg <- overline(routes, attrib = "n", fun = sum)
routes <- routes[order(routes$distance, decreasing = TRUE), ]
l_seg <- list()
# La fonction line_segment découpe un objet sf LINESTRING en segments réguliers
for(i in 1:nrow(routes)){
l_seg[[i]] <- line_segment(
routes[i, ],
segment_length = max(routes$distance) / 20)
}
pal <- mf_get_pal(nrow(l_seg[[1]]), rev = TRUE, palette = "Blues 3")
om_map(res, title = "Routes", theme = "dark")
for(i in nrow(routes):1){
mf_map(l_seg[[i]], col = pal, lwd = 1.6, add = TRUE)
}
leg(type = "cont", val = c(min(routes$distance), 1000, max(routes$distance)),
pal = pal, pos = "left", title = "Distance (m)")Le package osrm propose aussi des isochrones, à l’aide de la fonction osrmIsochrone().
o <- data.frame(X = st_coordinates(pt)[, 1],
Y = st_coordinates(pt)[, 2])
isos <- osrmIsochrone(loc = o,
breaks = seq(0,30,5),
res = 60,
osrm.profile = "foot")
saveRDS(isos, "data/isos.rds")isos <- readRDS("data/isos.rds")
isos <- st_transform(isos, crs = "EPSG:3857")
isos <- st_intersection(isos, res$zone)
om_map(res, title = "Isochrones autour de la mairie de Compiègne")
mf_map(isos, type = "typo", var = "isomax", alpha = .6, pal = "Geyser",
border = "white", leg_pos = NA, add = TRUE)Bonus : la tournée des pizzeria
Je me trouve à la mairie de Compiègne et souhaite goûter toutes les pizzas de la ville en un minimum de temps. Avec la fonction osrmTrip(), c’est possible !
Cette fonction répond au problème du voyageur de commerce, qui consiste à déterminer le plus court circuit passant par chaque point une seule fois.
# Coordonnées des pizzeria
pizz <- poi[!is.na(poi$cuisine) & poi$cuisine == "pizza" , "geometry"]
# Ajouter le point de départ : la mairie
pizz <- rbind(pt[1, "geometry"], pizz)
pizz$id <- 1:nrow(pizz)
trip <- osrmTrip(loc = pizz,
overview = "full",
osrm.profile = "foot")
saveRDS(trip, "data/trip.rds")trip <- readRDS("data/trip.rds")
trips <- trip[[1]]$trip |>
st_transform(crs = "EPSG:3857")
# Pour les labels : récupérer le début de la ligne
start <- lwgeom::st_startpoint(trips)
start <- do.call(rbind, start)
start <- data.frame(id = rownames(trips), X = start[,1], Y = start[,2])
start <- st_as_sf(start, coords = c("X", "Y"), crs = "EPSG:3857")
pizz <- poi[!is.na(poi$cuisine) & poi$cuisine == "pizza" ,]|>
st_transform(crs = "EPSG:3857")
om_map(res, title = "La route de la pizza")
mf_map(res$zone, col = "#a83c0a33", alpha = .1, border = NA, add = TRUE)
mf_map(res$zone, col = NA, border = "#d7b578", lwd = 15, add = TRUE)
mf_map(trips, col = "black", lwd = 2, add = TRUE)
mf_map(start[start$id != 1, ], pch = 20, cex = 2, col = "#496d0e", add = TRUE)
mf_map(pt_reproj, pch = 24, cex = 1.3, col = "darkblue", lwd = 2, add = TRUE)
mf_label(start, var = "id", cex = 1.1, overlap = TRUE, pos = 4, halo = TRUE)